# fitting-testes.R - testes de hipóteses

# Exemplo de uso:
#
# source('testes-pvalue.R')
# source('dados.R')
# Fitting(Danos(3))

library(actuar)
library(nortest)
library(MASS)

## Testes de hipóteses

TesteNormalShapiro <- function(dados)
	data.frame(pvalue = shapiro.test(dados)$p.value)

TesteNormalLillie <- function(dados)
	data.frame(pvalue = lillie.test(dados)$p.value)

TesteGamma <- function(dados) {
	E <- mean(dados) ; V <- var(dados)
	alpha <- (E^2)/V
	beta <- E/V
	if(alpha > 0 && beta > 0)
		return(data.frame(
			pvalue = ks.test(dados, y="pgamma", alpha, beta)$p.value,
			vals = paste(
				'alpha:', format(alpha, digits=6),
				'beta:', format(beta, digits=6))))
	0
}

TesteChisq <- function(dados) {
	E <- mean(dados)
	data.frame(
		pvalue = ks.test(dados, y="pgamma", E/2, 0.5)$p.value,
		vals = paste(
			'k:', format(E, digits=6)))
}

# caso particular de gamma, em q alpha=1

TesteExp <- function(dados)
{
	beta <- mean(dados)/var(dados)
	if(beta > 0)
		return(data.frame(pvalue = ks.test(dados, y="pgamma", 1, beta)$p.value))
	0
}

TesteUniforme <- function(dados)
{
	binomio <- 4 - (4*(mean(dados)^2) - 12*var(dados))
	if(binomio >= 0)
		for(i in c(-1,+1)) {
			c <- 2 + binomio*i
			if(c >= 0)
				for(j in c(-1,+1)) {
					b <- sqrt(c)*j
					a <- 2*E - b
					if(a < b)
						return(data.frame(
							pvalue = ks.test(dados, y="punif", a, b)$p.value,
							vals = paste(
								'a:', format(a, digits=6),
								'b:', format(b, digits=6))))
				}
		}
	0
}

TestePareto <- function(dados)
{
	E <- mean(dados) ; V <- var(dados)
	c <- (E^2) / V
	binom <- 4 - 4*1*(-c)   # binomio discrimante
	if(binom >= 0)
		for(i in c(-1,+1)) {
			alpha <- (2 + i*sqrt(binom)) / 2
			if(alpha > 0) {
				x0 <- E * (alpha-1) / alpha
				if(x0 > 0)
					return(data.frame(
						pvalue = ks.test(dados, y="ppareto", alpha, x0)$p.value,
						vals = paste(
							'alpha:', format(alpha, digits=6),
							'x0:', format(x0, digits=6))))
			}
		}
	0
}

## Batch de testes

TestesFn <- c(
	TesteNormalShapiro, TesteNormalLillie,
	TesteGamma, TesteExp, TesteChisq, TestePareto)
TestesNome <- c(
	"Normalidade Shapiro-Wilk", "Normalidade Lilliefors",
	"Gamma", "Exponencial", "Chi-square", "Pareto")

BatchTestes <- function(dados, descr)
{
	i <- 1
	for(fn in TestesFn) {
		res <- fn(dados)
print(TestesNome[i])
print(res)
		if(res$pvalue >= 0.05) {
			out <- paste(descr, 'com', TestesNome[i],
				', p-value:', format(res$pvalue, digits=6))
			if(length(res$vals) > 0)
				out <- paste(out, '-', res$vals)
			print(out)
		}
		i <- i+1
	}
}

Fitting <- function(Dados, sqmax = 0)
{
	print("nao rejeitado por:")
	BatchTestes(Dados, 'Dados')
	BatchTestes(log(Dados), 'log(Dados)')
	# NOTA: não vale a pena calcular log(Danos) com outras bases
	BatchTestes(1/Dados, "Dados^-1")
	if(sqmax > 1)
		for(b in c(2:sqmax))
			BatchTestes(Dados^(1/b), paste("raiz(Dados) base ", b))
}

